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Techniques are developed for projecting the solutions of symmetric hyperbolic evolution systems 
onto the constraint submanifold (the constraint-satisfying subset of the dynamical field space). 
These optimal projections map a field configuration to the "nearest" configuration in the constraint 
submanifold, where distances between configurations are measured with the natural metric on the 
space of dynamical fields. The construction and use of these projections is illustrated for a new rep- 
resentation of the scalar field equation that exhibits both bulk and boundary generated constraint 
violations. Numerical simulations on a black-hole background show that bulk constraint violations 
cannot be controlled by constraint-preserving boundary conditions alone, but are effectively con- 
trolled by constraint projection. Simulations also show that constraint violations entering through 
boundaries cannot be controlled by constraint projection alone, but are controlled by constraint- 
preserving boundary conditions. Numerical solutions to the pathological scalar field system are 
shown to converge to solutions of a standard representation of the scalar field equation when con- 
straint projection and constraint-preserving boundary conditions are used together. 



I. INTRODUCTION 

The exponential growth of constraint violations in the 
evolutions of black-hole spacetimes is probably the most 
critical problem facing the numerical relativity commu- 
nity today. The evolution equations of any self-consistent 
evolution system with constraints (including Einstein's) 
ensure that if the constraints are satisfied identically on 
an initial spacelike surface, they will remain satisfied 
within the domain of dependence of that surface. This 
does not mean that small initial violations of the con- 
straints will remain small, or that constraint violations 
will not flow into the computational domain through 
timelike boundaries. On the contrary, experience has 
shown that constraint violations seeded by roundoff or 
truncation level errors in the initial data tend to grow 
exponentially in the numerical evolutions of black-hole 
spacetimes (see e.g., 0,0,01). At present these constraint 
violating instabilities are the limiting factor preventing 
these numerical simulations from running for the desired 
length of time. Finding ways to control the growth of 
these constraints is therefore our most urgent priority. 

Recent work has demonstrated numerically that con- 
straint violations that flow into the computational do- 
main through timelike boundaries can be controlled ef- 
fectively by the use of special constraint preserving 
boundary conditions Q H H Q, ||. A number of 
groups have constructed such boundary conditions for 
various representations of the Einstein evolution sys- 
tem |!E0l![lllp[ll[30HQl. However, 
constraint violations in many evolution systems (includ- 
ing Einstein's) are driven by bulk terms in addition to 
boundary terms in the constraint evolution equations. In 
this paper we demonstrate that such bulk generated con- 
straint violations cannot be controlled effectively through 
the use of boundary conditions alone. Alternative meth- 



ods of controlling the growth of constraints are still re- 
quired in such systems. 

The most widely used method of controlling the growth 
of constraints in the Einstein evolution system is called 
fully constrained evolution. In this method, which is of- 
ten applied to spherical or axisymmetric problems, sym- 
metry considerations are used to separate the dynamical 
fields into those that are determined by solving evolution 
equations and those that are determined by enforcing the 
constraints at each time step [ill 111 HI IM l2ll I2al23| . 
In 3D problems without symmetry there is no obvious 
way to perform such a separation in a general coordinate 
system; however, fully constrained 3D methods based on 
spherical coordinates have yielded promising results [24| . 
Various groups have studied a closely related method, 
constraint projection, which can be used for general 3D 
evolutions in any coordinate system. The idea is to use 
the evolution system to advance all of the dynamical 
fields in time, and then at each time step (or whenever 
the constraints become too large) to force the solution 
back onto the constraint submanifold by solving the con- 
straint equations (for the conformal factor and the longi- 
tudinal part of the extrinsic curvature in the case of the 
Einstein system). The first preliminary results obtained 
with this constraint projection technique have been mod- 
erately successful [24], |25|, |26| . Constraint projection has 
not gained widespread use in 3D simulations, however, 
due in part to the traditionally high cost of solving the 
elliptic constraint equations. Difficult questions also re- 
main unanswered about the proper boundary conditions 
to impose on the constraint equations, for example at 
black hole excision boundaries. Moreover, little attention 
has been given to the question of whether these projec- 
tions correctly map a field configuration onto (or near) 
the correct point of the constraint submanifold, i.e. the 
point through which the exact evolution of the system 
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would pass at that time. In particular, it is not clear 
whether the overall time evolution scheme — including the 
projections — remains consistent, stable, and convergent. 

The need to enforce constraints is a common feature of 
many problems in mathematical physics besides numer- 
ical relativity, and for many problems successful tech- 
niques have been developed to ensure that numerical so- 
lutions satisfy the needed constraints. Under mild as- 
sumptions on the constraints, the subset of the field 
space satisfying the constraint equations defines a for- 
mal differentiable manifold (a classical result due to 
Ljusternik |27}l. and the evolution of a dynamical sys- 
tem of ordinary (ODE) or partial differential equations 
(PDE) subject to constraints may be viewed as evolution 
on this submanifold. Constraint control methods for such 
systems are generally based on ideas from variational me- 
chanics, where the Lagrangian (whose stationary points 
describe the physical states of the system) is augmented 
with a sum of terms consisting of products of Lagrange 
multipliers and the constraints jH, US HO] • A necessary 
condition for a configuration point to be a solution of 
both the field equations and the constraint equations is 
that the augmented Lagrangian be stationary with re- 
spect to variations in both the fundamental fields and 
the Lagrange multipliers jS^. The additional terms in 
the augmented Lagrangian involving Lagrange multipli- 
ers can be viewed as forcing the dynamics to remain on 
the constraint submanifold. 

These augmented Lagrangian techniques are the basis 
of well-studied numerical methods for controlling con- 
straint violations in ODE systems. Many ODE systems 
are subject to algebraic constraints which must be pre- 
served as the solution evolves. For such systems there 
exist numerical integration techniques that enforce these 
algebraic constraints exactly, and that also conserve vari- 
ous important properties of the ODE solution (e.g. time- 
reversibility and symplcctic structure). These numerical 
techniques are derived by adding to the ODEs terms cho- 
sen to make a suitable augmented Lagrangian for the sys- 
tem stationary |32, IH, 03 . The resulting numer- 
ical schemes, referred to as "step-and-project" methods, 
can be thought of as standard time integration steps fol- 
lowed by projections. First a preliminary step is taken 
forward in time using a standard numerical scheme, af- 
ter which the solution will generally not satisfy the con- 
straint equations. Then the solution from the prelimi- 
nary step is corrected using a formal (optimal, or near- 
est point) projection back onto the constraint submani- 
fold. This projection step typically involves solving alge- 
braic equations. Unlike the simple constraint projection 
methods used so far in numerical relativity, "step-and- 
project" numerical methods for constrained systems are 
well studied and well understood. It has been shown 
that they retain the consistency and stability proper- 
ties of the original one-step method on which they are 
based, and they generally have the same convergence 
properties |35| . These techniques are immediately appli- 
cable to constrained PDE systems that are discretized 



in space to produce constrained ODE systems (as we 
do, see Sec. llVf) : and numerical methods based on aug- 
mented Lagrangians for PDE systems have also been de- 
veloped [39.137^ 

In this paper we apply these augmented variational 
techniques to obtain equations that project solutions of 
constrained hyperbolic evolution systems onto the con- 
straint submanifold of the appropriate dynamical field 
space. We construct projections that are optimal, in 
the sense that they map a given field configuration to 
the "nearest" point on the constraint submanifold. We 
use the natural metric, the symmetrizer, that exists in 
any symmetric hyperbolic evolution system to define dis- 
tances on the space of fields. Hence this optimal projec- 
tion is the one that minimizes this symmetrizer distance 
(typically called the energy) between a given field con- 
figuration and its projection. The general formalism for 
constructing such optimal projections for constrained hy- 
perbolic evolution systems is described in Sec. ITT1 

We illustrate these optimal constraint projection ideas 
in Sec. IHII bv deriving the optimal projection for a new 
symmetric-hyperbolic representation of the scalar field 
equation on a fixed background spacetime. This scalar 
field system has the interesting property that it suffers 
from constraint violations driven both by bulk terms as 
well as boundary flux terms in the equations. (And so 
this system serves as a good model of the pathologies 
present in the Einstein system.) The optimal projection 
for this scalar field system is determined by solving a 
certain elliptic PDE. In Sec. IIVI we test these optimal 
projection techniques by studying numerical solutions to 
this scalar field system on a fixed black-hole background 
spacetime. In particular we demonstrate that constraint 
preserving boundary conditions are necessary, but not 
sufficient, to control the growth of constraints in this 
pathological scalar field system. We demonstrate that 
constraint projection succeeds in producing convergent 
constraint-satisfying solutions, but only if constraint pre- 
serving boundary conditions are used as well. These tests 
also illustrate that the projections are best performed 
at fixed time intervals (AT w 2M for this problem) 
rather than after each time step. And we show that 
the computational cost of solving the constraint projec- 
tion equations for this system (using our spectral ellip- 
tic solver 38] ) is a very small fraction (below 1% for 
the resolution needed to achieve roundoff level accuracy) 
of the total computational cost of evolving this system. 
The symmetrizer metric for this model scalar field system 
(like many hyperbolic evolution systems) is not unique; 
so the projections defined in terms of the symmetrizer are 
not unique. Nevertheless, we demonstrate for the model 
scalar field system that numerical evolutions based on 
these different projections all converge to the same so- 
lution. The rate of this convergence is not the same for 
all projections, however, and we find an "optimal" pro- 
jection for this system that maximizes this convergence 
rate. 
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II. OPTIMAL CONSTRAINT PROJECTION 

Our objective is to construct a projection operator that 
maps a given field configuration to the nearest constraint- 
satisfying configuration (the nearest point on the con- 
straint submanifold) . That is, we wish to map an initial 
point u a in the field configuration space to a new point 
u a that satisfies a set of constraint equations: 

c A {u a ) = 0. (1) 

(We use Greek indices to label individual components of 
the dynamical fields, and upper case Latin indices to la- 
bel the individual components of the constraints.) To 
find the optimal projection we also need to have a dis- 
tance measure between field points. We define the needed 
measure in terms of a symmetric positive-definite metric, 
S a p, on the dynamical field space. The distance between 
field points is then defined as 

\\Su\\ 2 = J S al3 {u a ~u a ){u p -u^x. (2) 

Building on the augmented variational techniques com- 
monly used to construct step-and-project constraint con- 
trol schemes in other areas of numerical analysis |3lL l33l 
l35l |. we are now prepared to construct the optimal pro- 
jection map. We introduce a Lagrangian density C that 
consists of the distance between the given field config- 
uration u a and its projection u a , plus the products of 
the constraints with Lagrange multipliers. Thus we in- 
troduce the Lagrangian density, 

C = S aP (u a -u a )(u (3 -u f) ) + \ A c A . (3) 

The stationarity of the Lagrangian (the volume integral 
of this Lagrangian density) with respect to variations 
of the Lagrange multipliers A^ enforces the constraints, 
while stationarity with respect to variations of the fields 
u a are necessary conditions for the projection to mini- 
mize the distance to the constraint submanifold. 

The optimal projection procedure outlined above could 
be carried out using any definition of the distance be- 
tween field points, e.g. using any positive definite metric 
S a p on the space of fields. For a particular problem this 
distance measure should be chosen to be the natural mea- 
sure associated with that problem. Our primary interest 
here is the construction of projections for constrained hy- 
perbolic evolution systems. So we will focus our attention 
on fields u a that satisfy a first-order evolution equation 
of the form 

d t u a + A koi ^d h vP = F a . (4) 

We use lower case Latin indices to label spatial coor- 
dinates x k , d t = d/dt to denote time derivatives, and 
dk = d/dx k to denote spatial derivatives. Such systems 
are called symmetric hyperbolic if they have a positive 
definite metric S a p on the space of fields (typically called 



the symmetrizer) that symmetrizes the characteristic ma- 
trices: 

S ai A kl p = A a p — Ap a . (5) 

The well-posedness of the initial value problem for linear 
symmetric-hyperbolic evolution systems is demonstrated 
by establishing bounds on the square-integral norm of 
the dynamical fields defined with this symmetrizer met- 
ric |39l |40| . This metric defines the meaningful measure 
on the dynamical field space for symmetric-hyperbolic 
systems, so this is the appropriate measure to use for 
constructing optimal constraint projections for these sys- 
tems. Most hyperbolic evolution systems of interest in 
mathematical physics (including many representations of 
the Einstein system) are symmetric hyperbolic, and so we 
limit our consideration here to systems of this type. 

In Sec. lIIII we use the procedure outlined above to con- 
struct explicitly the optimal projection for the relatively 
simple case of the scalar wave equation on a curved back- 
ground spacetime. But before we focus on that special 
case, we take a few (rather more abstract) steps in the 
construction of this projection for the general case. To 
do this we assume that the constraints c A are linear in 
the derivatives of the dynamical fields: 

c A = K kA (i d k v} 3 + L A , (6) 

where K kA p and L A may depend on u a but not its 
derivatives. The constraints have this general form in 
many evolution systems of interest {e.g., the Einstein 
system, the Maxwell system, the incompressible fluid sys- 
tem). In this case we can explicitly compute the varia- 
tions of the Lagrangian density defined in Eq. (J3J : 

^6u a =5u a {2S a p(uP -uP)-d k (\ A K kA a ) 
+\ A (d Q K kA pd k uP + d a L A )} 

+d k (X A K kA a 6u a ), (7) 
^6X A = c A 5X A . (8) 

Here we use the notation d a = d/du a to denote deriva- 
tives with respect to the fields. We have also assumed 
that the symmetrizer S a p may depend on u a but not u a . 
We wish to find the stationary points of this Lagrangian 
with respect to arbitrary variations in the fields u a and 
the Lagrange multipliers A^. Stationarity with respect 
to the variations of these quantities (that vanish on the 
boundaries) implies that 

= u a ~ u a ~ ^S ap d k {X A K kA p) 

+^X A S^{dpK kA 1 d k v? + d L A ), (9) 
= c A = K kA p d k u + L A (10) 

at each interior point, and stationarity with respect to 
boundary variations implies that 

= n k \ A K kA p (11) 
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at each boundary point, where n k is the outward directed 
unit normal to the surface. We use the notation S a/3 to 
denote the inverse of S a p. The general idea is to use 
Eqs. ® and (|10l) . with appropriate boundary conditions 
(such as those provided by Eq. m]), to determine the 
field configuration u@ and the Lagrange multipliers Xa 
for any given field point u a . If u a and Xa satisfying these 
equations can be found, then we are guaranteed that the 
field u a is the constraint-satisfying solution nearest the 
point u a as desired. We do not know whether these equa- 
tions always admit solutions in the general case. So in 
Sec. IIIII wc study in detail this optimal projection for 
the simple case of the scalar field equations on a fixed 
background spacetime. We show that solutions to the 
optimal projection equations always exist and are rela- 
tively easy to compute numerically in this simple case. 
And in Sec. IIVI we show that this optimal projection is 
very effective in controlling the growth of constraints for 
the scalar field system. 



III. SCALAR FIELDS IN CURVED SPACETIME 

In this section we examine in some detail the scalar 
wave system on a fixed background spacetime. In 
Sec. IIII Al we review the standard treatment of this sys- 
tem, and then modify it so that it exhibits bulk generated 
constraint violations in addition to the boundary gener- 
ated violations present in the standard system. This new, 
more pathological, symmetric-hyperbolic scalar field sys- 
tem now serves as a good model of the constraint violat- 
ing problems inherent in the Einstein system. We con- 
struct constraint preserving boundary conditions for this 
system in Sec. IIII Bl and the optimal projection map for 
this system in Sec. IIII d following the procedure outlined 
in Sec. El 



A. Modified Scalar Wave System 

The standard scalar field equation on a fixed back- 
ground spacetime is 

VV^V - 0, (12) 

where ip represents the scalar field and V M the covari- 
ant derivative associated with the background spacetime 
metric. We represent the background spacetime metric 
in terms of the usual 3 + 1 splitting: 

ds 2 = -N 2 dt 2 + giJ (dx l + N l dt){dx J + N 3 dt), (13) 

where the lapse N and the spatial metric gij are assumed 
to be positive definite, while the shift N 1 is arbitrary. 
The equation for the scalar field ip, Eq. (|12fl . can be re- 
expressed as a first-order evolution system in the stan- 



dard way (see e.g. Ref. [ill]'): 

d t t(j - N k d k tp = -ATI, (14) 

<9 t IT - N k d k Il + Ng ki d k $i = NJ t <fr l + NKU, (15) 
d t ^ l - N k d k <Pi + Nd t U = -UdiN + QjdiN 1 . (16) 

The field <I>i represents the spatial gradient diip, and II 
represents the time derivative of ip (and is defined pre- 
cisely by Eq. ^1]). The auxiliary quantities K (the trace 
of the extrinsic curvature) and J 1 depend only on the 
background spacetime geometry, and are defined by 

J 1 = -N^g-WdjiNgWgV), (17) 

K = -A^V^W 72 -^ 172 ^')]- ( 18 ) 

Solutions to the first-order system, Eqs. (|14l> -l|16 [l . are 
also solutions to Eq. I|12fl only if the constraints are sat- 
isfied: = c A = {C l ,C i] }, where 

C t = diip-Qi, (19) 

Cy = (20) 

Although the second constraint, = 0, follows from 
the first, Ci = 0, the converse is not true. Hence we 
include both constraints in the analysis here. Note that 
both constraints are necessary to construct a first-order 
hyperbolic evolution system for the constraint quantities 
(discussed below, Eqs. (|2*5|) and P0)l). Note also that the 
analogues of both constraints play essential roles in first- 
order hyperbolic formulations of Einstein's equations. 

We now generalize the evolution system, Eqs. I|14|) - 
(|16[) . somewhat by adding multiples of the constraint Ci 
to Eqs. O and (JTSJ: 

d t il> - N k d k i> = -NU + ^Ck, (21) 

d t <$>i - N k d k <Pi + NdiU = -UdiN + <^ j d l N j 

+j 2 Nd, (22) 

where 71 and 72 are arbitrary constants. The constraint- 
satisfying solutions to these equations are the same as 
those of the original system; but as we shall see, the con- 
straint violating properties of the new system are signif- 
icantly different from those of the original. Substituting 
the definition of Ci in Eqs. 1121 [I and l|22l) gives us new 
evolution equations for ip and 

d t ^-(l + ll )N k d k ^ = -NU- ll N k ^ k , (23) 
d t ^ - N k d k <S>i + NdiU - 72 = 

- UdiN + ^> 3 d l N j - 7 2 A^. (24) 

The first-order system that represents the scalar wave 
equation, Eqs. I|15|) . H23|) . and (|24|l . has the standard first 
order form, 

3 t u a + A ka p d k u p = F a , (25) 

where u a — {tp, II, <i>i}. Systems of this type are 
called symmetric hyperbolic if there exists a symmetric 
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positive-definite tensor S a p on the space of fields that 
symmetrizes the characteristic matrices A ka p\ 

S ai A kl /3 = A k afJ = Ap a . (26) 

The most general symmetrizer for our new scalar wave 
system is (up to an overall factor), 

ds 2 = S a pdu a du^ , 

= A 2 d^ 2 - 2 l2 dij dU + dU 2 + gVdQidQj, (27) 

where A is an arbitrary non- vanishing function. This S a p 
symmetrizes the characteristic matrices A a/3 so long as 
7i72 = 0. Thus we must take at least one of these pa- 
rameters to be zero for our new system to be symmetric 
hyperbolic. This symmetrizer is positive definite when- 
ever 

A 2 > 7 |. (28) 

In this case S a p provides a dynamically meaningful mea- 
sure of the distance between field configurations, which 
we use to define our optimal constraint projection oper- 
ator in Sec. cnn 

The evolution of the constraints follows from the prin- 
cipal evolution system, Eqs. lfT5|l . J2HJ), and (f2"H> : 

dtCi-il + ^Ctfd =2rf X N!C ii -t l NCi, (29) 
d t Cij - L n Cij = -TaJVCy - 72%^ N, (30) 

where represents the Lie derivative along the shift 
vector N l . If the constraints are satisfied at some ini- 
tial time, then these equations guarantee (at least at 
the analytical level) that the constraints remain satis- 
fied in the domain of dependence of the initial data. 
These equations also show that any constraint violations 
in this system will be advected along a congruence of 
timclike curves. Constraint violations can therefore flow 
into the computational domain if these curves intersect 
the boundaries. And like the Einstein evolution sys- 
tem, these equations also contain bulk terms that am- 
plify any existing constraint violations. When 71 = 
we see that Eq. <|29(1 implies that the constraint Ci has 
the simple time dependence Ci(r) = Ci(0)e _72T , where r 
measures proper time as seen by a hypersurface orthog- 
onal observer. Whenever 72 < this constraint grows 
exponentially, and in this case the modified scalar wave 
system serves as a good model of the constraint viola- 
tions in the Einstein system. (Constraint violations of 
all wavelengths grow exponentially in this system, and so 
it may be even more pathological than the Einstein sys- 
tem where constraint violating instabilities are typically 
dominated by long wavelength modes @j|3].) Conversely, 
if 72 > then this modified scalar wave system exponen- 
tially suppresses any residual constraint violations that 
may be present in the initial data. This latter prop- 
erty suggests that analogous terms could be introduced 
to control some of the bulk constraint violating terms in 
the Einstein system. 



B. Constraint Preserving Boundary Conditions 

Boundary conditions for hyperbolic evolution systems 
are defined in terms of the characteristic fields of these 
systems, so we must construct these fields for our modi- 
fied scalar wave system. The characteristic fields are de- 
fined with respect to a spatial direction at each point, rep- 
resented here by a unit normal co- vector field n k . For the 
purposes of imposing boundary conditions, the appropri- 
ate n k is the outward-pointing normal to the boundary. 
Given a direction field n k we define the left eigenvectors 
e a a of the characteristic matrix rikA ka p by 

e a a n k A ka p = v^ a )e a p, (31) 

where V( a ) denotes the eigenvalue (also called the charac- 
teristic speed). The index a labels the various eigenvec- 
tors and eigenvalues, and there is no summation over a in 
Eq. jSJ . Since we are interested in hyperbolic evolution 
systems, the space of eigenvectors has the same dimen- 
sion as the space of dynamical fields, and the matrix e a p 
is invertible. The projections of the dynamical fields u a 
onto these left eigenvectors are called the characteristic 
fields u a : 

u & = e & pu (i . (32) 

At each boundary point, boundary conditions must be 
imposed on any characteristic field having negative char- 
acteristic speed, < 0, at that point ^2, E3- We 
refer to fields with V/ty < as the incoming characteris- 
tic fields at that point. Conversely, those characteristic 
fields having non-negative characteristic speeds (the out- 
going fields) must not have boundary conditions imposed 
on them there. 

The characteristic fields for the symmetric hyperbolic 
representations (7172 = 0) of the scalar wave system are 
the quantities u a = {Z 1 , Z 2 , U 1 ^}: 

Z 1 = V, (33) 
Z 2 = P k ^ k , (34) 
U 1 * = n ± n fe $ fc - 72^, (35) 

where P k i — 5 k i — n k rii, n k = g ■'rij, and n k n k = 1. 
The fundamental fields u a can be reconstructed from the 
characteristic fields u a by inverting Eq. (|32|) : 

tl> = Z 1 , (36) 

n = ±(u 1+ + u 1 -) + l2 z\ (37) 
*< = h(u 1+ -U^m + Zf. (38) 

The characteristic field Z 1 propagates with speed 
—"/in k N k /N, the field Z 2 with speed 0, and the fields 
C/ 1± with speeds ±1 relative to the hypersurface orthog- 
onal observers. The coordinate characteristic speeds of 
these fields are — (1 + ji)n k N k , —n k N k and —n k N k ±N 
respectively. 

At each boundary point, boundary conditions are re- 
quired on each characteristic field whose coordinate char- 
acteristic speed is negative at that point. The field C/ 1_ , 
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in particular, requires a boundary condition on all time- 
like boundaries. For the standard representation of the 
scalar field system, Eqs. lfT^jl ~ (|T?))l . the boundary condi- 
tion U 1 ^ = Id — n k $k = is used to ensure (approx- 
imately) that no scalar waves enter the computational 
domain. We wish to enforce this condition on our gen- 
eralized scalar field system, Eqs. <|15fl . 1221 j an d (1241) . 
in such a way that the physical (constraint satisfying) 
solutions are the same for all values of the parame- 
ters 7! and 7 2 . Since U 1 ^ depends on 72, Eq. (|35fl . 
the proper boundary condition must also depend on 72 : 
U 1 ^ + 72V' = n — n fe $fc = 0. Thus the appropriate 
boundary condition to impose on U 1 ^ is U 1 ^ = —72^- 
The freezing form of this boundary condition (as used in 
our code) is, 



the constraint submanifold. Following the procedure out- 
lined in Sec. ^ we construct a Lagrangian density, 



dtU 1 = ~j 2 d t tp- 



(39) 



For boundary conditions on the fields Z 1 and Zf (when 
necessary), we explore two choices: One is the freezing 
boundary condition d t Z l — d t Z 2 — 0. In Sec.|^we show 
that this boundary condition allows constraint violations 
to enter the computational domain through the bound- 
aries. Therefore, we also explore conditions that prevent 
this influx of constraint violations: When the fields Z 1 
and/or Z 2 require boundary conditions, we set 



dtZ 1 = 7V fc $ fe -ivn, 
d t zf = P k id k d t i>. 



(40) 
(41) 



Equation (|40|l is based on Eq. l|14fl combined with 
Eq. (|19|l . while Eq. i|41|) is derived from the time- 
derivative of Eq. I|19|l . We note that with the choice 
71 = —1, the field Z 1 never requires a boundary condi- 
tion. We also note that the term dttp that appears on 
right side of Eqs. l|3*9"|) and i|4*T1l must be evaluated us- 
ing the appropriate expression for dt4> — dtZ 1 on this 
boundary: Eq. (|4(J|) when Z 1 requires a boundary condi- 
tion, or Eq. 12M[) when no boundary condition is required. 
In Sec. II VI we compare numerically the results of using 
these constraint preserving boundary conditions with the 
use of the freezing boundary conditions dtZ 1 = dtZ 2 = 
on these fields. 



C. Optimal Constraint Projection 

The idea is to use the full evolution system, Eqs. (1151) . 
Q23JI. and l|24|l . to evolve initial data forward in time 
an amount AT and then (when the constraint viola- 
tions become too large) to project this solution back 
onto the constraint submanifold in some optimal way. 
Let u a = {-0,11, $j} denote the solution obtained di- 
rectly from this free evolution step. This solution u a 
may not satisfy the constraints because roundoff or trun- 
cation level constraint violations have been amplified, or 
constraint violations have flowed through the boundaries. 
Thus we wish to project u a in an optimal way back onto 



A 2 (v>-^) 2 - 2 72 (v-^)(ii-n) 
-(n - n) 2 + - $;)($,- - $, ) 



= 9 1/2 



(42) 



using the symmetrizer S a p of the hyperbolic evolution 
system, Eq. (|27|l . and the Lagrange multipliers = 
{A\ A iJ }. The stationary points of the Lagrangian, 



L= Cd 



(43) 



with respect to variations in u a and A ,4 define the opti- 
mally projected field configuration u a . We have included 
the multiplicative factor g 1 / 2 = (det gij) 1 ^ 2 in Eq. l|4^|) 
to ensure that L is coordinate invariant. 

The scalar field constraint Lagrangian density, 
Eq. l|42(l . has the following variations: 

§5V =2 5 1 / 2 [A 2 (^-0)-72(n-n)]^ 

-d^X^ + W'XS^), (44) 

§<m =2 ff 1 / 2 [n-n-7 2 (0-0)] ( 5n, (45) 

j§-5^ = [2 5 1 /2^ ($ ._$. ) _ 5 i/2 A i 

-d^g 1 ' 2 ^)^, + 9 l (.g 1 / 2 A l ^$ J ),(46) 
faSX* =g 1 t 2 {d^-<S> l )5X\ (47) 
faSX* ;/' 2 ^,<MA". (48) 

We require that the Lagrangian L from Eq. (|43|l be sta- 
tionary with respect to all variations in the dynamical 
fields Su a — {Sip, <5II, <5$i} (including those that do not 
vanish on the boundaries) as well as all variations in the 
Lagrange multipliers SXa = {SX l ,SX lJ }. From Eqs. l|4^|) ~ 
(gSJ), it follows that 

V; = ^ + 72A- 2 (n-n) + iA- 2 5 - 1 /2^( ff i/2y ); (49) 
n = n + 7 2 (v>-V), (50) 

= Si + kQaX + ig- 1/2 g t] d k (9 1/2 ^ 3 ), (51) 

and Eqs. (|47|l and l|48() imply that the projected solution 
satisfies the constraints. We now solve Eq. 1|> for A 1 , 
substitute it into Eq. jHfll , and use Eqs. iflTjl and JHOJ, 

to obtain the following equation for ip, 

V 4 V^ - (A 2 - 72 2 )^ = V% - (A 2 - 7 2 2 )0, (52) 

where Vi represents the spatial covariant derivative that 
is compatible with g^j . In deriving this equation we have 
also used the fact that the term did^g 1 ^ 2 X kl ) vanishes 
identically because A y is antisymmetric. Equation l|52|) 
is just the covariant inhomogeneous Helmholtz equation. 
We note that the parameters must satisfy the condition 
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A 2 — 7| > f° r the evolution system to be symmet- 
ric hyperbolic. Solving Eq. 1)52)1 determines the optimal 
projection ip; the optimal II is determined from Eq. (jHOJ, 

n = fi + l2 {4,-i>)- (53) 

and the optimal $^ is obtained by enforcing the con- 
straint, 

$ 4 = drf. (54) 

We note that the Lagrange multiplier A iJ does not play 
any essential role in this analysis: we could just as well 
have set A 4 - 7 = and still obtained the same projection. 
This makes sense, because the constraint C{j is really a 
consequence of the constraint Ci in this case. 

The evolution equations for II and <&i, Eqs. (|15|l and 
1)16)1 . decouple from the larger scalar field evolution sys- 
tem, Eqs. l|T5)l. and when 72 = 0. It is some- 
times of interest to consider the properties of this smaller 
system, Eqs. (11511 and (|lf>|> . subject to the single con- 
straint, Eq. 1)20(1 . The optimal constraint projection for 
this reduced system consists of Eqs. 1)51))) and 1)51)1 (with 
A 4 = 72 = 0), together with the single constraint equa- 
tion d[i$>j] = 0. This constraint equation implies that 
&i = dii/j for some scalar function ip. Inserting this ex- 
pression for <I>i in Eq. 1)51)) . multiplying by g^^g^, and 
taking the divergence, we obtain the following equation 
for ip, 

V'V.V = V%. (55) 

In deriving this equation we have used the fact that the 
term did k (g 1 ^ 2 \ kl ) vanishes identically because A 4 - 7 is an- 
tisymmetric. The optimal projection in this reduced sys- 
tem then sets II = II and <I>i = diip, where ip is the 
solution to Eq. l(55*|) . We note that Eq. lf53)> is just the 
A 2 — 7I = limit of the original projection Eq. 1)52)) . 

Unfortunately the optimal constraint projection for the 
scalar field system is not unique, because the parameter 
A in the symmetrizer metric is not unique. We have 
seen that taking the limit A 2 — >■ 7! is equivalent to ig- 
noring the evolution of the scalar field ip in constructing 
the optimal projection. Alternatively, the limit A — > 00 
corresponds to the simple projection ip = ip, II = II, 
and $i = diip. In this limit, no elliptic equation has to 
be solved, and the evolution of the field $i is effectively 
ignored when constructing the projection. We expect 
that the optimal choice of A will be one for which 1 / A 
corresponds to some characteristic length or time scale 
associated with the particular problem. We explore in 
Sec. UVCl the properties of these projection operators for 
a range of A, and show that an optimal value does exist. 
When 72 7^ the optimal choice seems to be A 2 = 27I, 
where I/I72I is the time scale on which the constraints 
are amplified. 

Finally, we must consider the boundary conditions for 
the projection equations that determine ip, i.e. Eq. 1)52)) 



or (|55|) . In general, boundary conditions for the pro- 
jection equations must satisfy two criteria: First, they 
must be consistent with boundary conditions imposed 
on the evolution equations, and second, the projection 
equations plus boundary conditions must not modify so- 
lutions that already satisfy the constraints. Typically, we 
enforce approximate outgoing wave boundary conditions 
on the evolution equations. For the case of the scalar 
wave equation, the approximate outgoing wave bound- 
ary condition, Eq. (J2HJ), sets U 1 ' = —72^ or equivalently 
n k §k = II on the boundaries (where n k is the outward 
directed unit normal). Since $^ = diip in these projected 
solutions, the appropriate boundary condition to impose 
on ip in Eq. 1)52)1 or 1)55)1 in this case would be 

n k d k iP = n = IL + l2 (iP~iP). (56) 

Alternatively we can derive boundary conditions for 
ip from the requirement that the boundary variations of 
the Lagrangian vanish. The divergence terms in Eqs. 1)44)) 
and 1)46)1 imply that 

= n k X k = n k \ kl , (57) 

on the boundaries for the scalar field system. A 
short calculation (using the fact that n k is proportional 
to a gradient, and X kl is antisymmetric) shows that 
nid k (g 1 / 2 X kl ) = 0, so we see from Eq. (|5l"1l that the nat- 
ural boundary condition is 

n k d k iP = n k $ k . (58) 

If the approximate outgoing wave boundary condition, 
n k <&k = ft, was used in the free evolution step, then the 
natural boundary condition Eq. 1)58)1 differs from Eq. 1561) 
by the term 72 (ip — ip). For the constraint projections 
described in Section llVI we impose the Robin boundary 
condition Eq. 1)561) on the solutions of Eq. ()52)) at the 
boundaries where U x ~ requires a boundary condition in 
the evolution step, and Eq. ()58)l on the solutions at all 
other boundaries (e.g. inside an event horizon). We note 
that the discrepancy between the natural and the physi- 
cal outgoing boundary condition could be eliminated by 
adding an appropriate boundary term to the constraint 
projection Lagrangian. 

IV. NUMERICAL RESULTS 

We have studied the effectiveness of the optimal con- 
straint projection methods developed in Sees. HT1 and IIIII 
for the case of a scalar field propagating on a fixed black- 
hole spacetime. For these simulations we use the Kerr- 
Schild form of the Schwarzschild metric as our back- 
ground geometry: 

2 M 

ds 2 = -dt 2 + (dt + dr) 2 + dr 2 + r 2 dtf. (59) 

r 

We express all lengths and times associated with these 
simulations in units of the mass, M, of this black hole. 
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Our computational domain consists of a spherical shell 
extending from r min = 1.9M (just inside the black-hole 
event horizon) to r max = 11. 9M. For initial data we use a 
constraint satisfying Gaussian shaped pulse with dipolar 
angular structure, 

ijj = 0, (60) 

II = Y w {6,v)e- (r - rof,w \ (61) 
$ 4 = 0, (62) 

with ro = 5M and w = 1M. The value of II is about 
2 x 1CT 21 at the outer boundary of our computational 
domain, below the level of double precision roundoff er- 
ror. 

For the remainder of this section we describe briefly 
the numerical methods used to solve this problem. Then 
in Sec. lIV Al we describe three numerical simulations de- 
signed to explore the effects of boundary conditions on 
the evolution of the constraints in these solutions. In 
Sec. IIVBI we describe two additional numerical simula- 
tions that illustrate the effectiveness of constraint projec- 
tion in controlling the growth of constraints. And finally 
in Sec. IIV CI we explore ways to optimize the use of the 
constraint projection method and measure its computa- 
tional cost. 

All numerical computations presented here are per- 
formed using a pseudospectral collocation method. Our 
numerical methods are essentially the same as those we 
have applied to evolution problems with the Einstein sys- 
tem [HH H 0], with scalar fields 41], and with the 
Maxwell system yj. Given a system of partial differen- 
tial equations 

d t u a {^t) = .F a [u(x,t),fl! i u(x,t)], (63) 

where u a is a collection of dynamical fields, the solution 
m q (x, t) is expressed as a time-dependent linear combi- 
nation of N spatial basis functions 0fc(x): 

N-l 

uft(x,t)= (64) 

k=0 

We expand each scalar function (ip and II) and the Carte- 
sian components of each vector ($ x , $ y , and <J> Z ) in 
terms of the basis functions T n (p)Yi m (9, (p), where Yi m 
are spherical harmonics and T n (p) are Chebyshev poly- 
nomials with 

2r J'max ^min /nr-\ 

P= • (65) 

We use spherical harmonics with I < £ max = 5 and a 
varying number of Chebyshev polynomials with degrees 
N r < 81. Spatial derivatives are evaluated analytically 
using the known derivatives of the basis functions: 

N-l 

d lU %(x, t) = J2 «*(Wk(x)- (66) 

k=0 



Associated with the basis functions is a set of N c colloca- 
tion points x.j . Given spectral coefficients (t) , the func- 
tion values at the collocation points u a (x; , f ) are com- 
puted by Eq. I|64|l . Conversely, the spectral coefficients 
are obtained by the inverse transform 

Nc-l 

"fc(*)= J2 mu%(x i ,t)i> k (x i ), (67) 

i=0 

where Wi are weights specific to the choice of basis func- 
tions and collocation points; thus it is straightforward to 
transform between the spectral coefficients u%(t) and the 
function values at the collocation points Ujy(x,,£). The 
partial differential equation, Eq. H63[) ■ is now rewritten 
using Eqs. l|6^|l - l|6*7|l as a set of ordinary differential equa- 
tions for the function values at the collocation points, 

d t u%{^t)=gf[u N {^,t)l (68) 

where Qf depends on u%(xj 7 t) for all j. This sys- 
tem of ordinary differential equations, Eq. (JBSJ, is in- 
tegrated in time using a fourth-order Runge-Kutta al- 
gorithm. Boundary conditions are incorporated into the 
rig ht side of Eq. @El using the technique of Bj0rhus |45j. 
The time step is typically chosen to be about one fifth 
the distance between the closest collocation points, which 
ensures that the Courant condition is well satisfied. This 
small time step is needed to reduce the time discretiza- 
tion error to the same order of magnitude as the spatial 
discretization error at radial resolution N r = 61. 

Elliptic partial differential equations, Eq. I|52|) or Q55p. 
are solved using similar pseu dosp ectral collocation meth- 
ods. As detailed in Ref. [3q|. we consider a mixed 
real/spectral expansion of the desired solution 

! max I 

ip(p n ,6,ct>) = J2 12 W»(M)> (69) 

1=0 m=-l 

where p n (for n = 0, . . . , N r — 1) are the collocation points 
of the Chebyshev expansion in (rescaled) radius p. Given 
a set of coefficients ipimm we can evaluate the residual of 
the elliptic equation and the residual of the boundary 
conditions using expressions like Eq. (|66|l : the require- 
ment that each Yj m component (for I < l max ) of this 
residual vanishes at the radial collocation points results in 
a system of algebraic equations for the coefficients ipimn- 
For the problem considered here these algebraic equa- 
tions are linear, and with suitable preconditioning are 
solved using standard numerical methods like GMRES. 
The elliptic solver is described in detail in Ref. (38|. 

We use no filtering on the radial basis functions, but 
apply a rather complicated filtering rule for the angular 
functions. When evaluating the right side of Eq. 168fl . 
we set to zero the coefficients of the terms with I = ^ max 
in the expansions of the scalars, dtip and <9 t IL The vec- 
tor dt$>i is filtered by transforming its components to a 
vector spherical harmonic basis, setting to zero the coef- 
ficients of the terms with I = l max in this basis, and then 
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FIG. 1: Constraint violations for evolutions with 71 = 72 = 0, 
freezing boundary conditions, and no constraint projections. 
Plotted are radial resolutions N r — 21, 31, . . ., 61; all curves 
lie on top of each other. 



transforming back to Cartesian components. The result 
-0 of each elliptic solve and the projected II (cf. Eq. [53]) 
arc filtered similarly. The projected <I>i is computed via 
Eq. H54JI from the filtered ip. We find no angular insta- 
bility, such as the one reported in Ref. when we use 
this filtering method. And we find no significant change 
in our results for this problem by increasing the value of 
f max beyond the value f mas = 5. 



A. Testing Boundary Conditions 

In this section we describe the results of three nu- 
merical simulations that explore the effects of bound- 
ary conditions on the evolution of the constraints in the 
scalar field system. First we evolve the initial data in 
Eqs. l|t)UI) - (|62|l using the standard representation of the 
scalar field system (71 =72 = 0), and using the stan- 
dard freezing boundary conditions on the incoming fields. 
We use no constraint projection in this initial simula- 
tion. At the inner boundary of the computational do- 
main, r — r m j n = 1.9M, all of the fields are outgoing and 
so no boundary condition is needed there on any of the 
fields. At the outer boundary, r — r max = 11.9M, the 
fields Z 1 , Zf and U l ~ are all incoming since the shift 
points out of the computational domain there: nkN k — 
2M/r. So we impose the freezing boundary conditions 
= dtZ 1 = d t Zf = dtU 1 - on these fields. The results 
of this first numerical simulation are depicted in Figs. 
and[3 

Figure ^ illustrates the evolution of the constraints, 



FIG. 2: Convergence plot for the evolution presented in Fig.Q 
Plotted are differences from the solution with radial resolution 
iV r = 81. 

which we measure using the quantity ||C(t)||, 

I |C(*)| 1 2 = / {C i C i + C ii &*)c?' i d 3 x i (70) 

divided by a suitable normalization. The constraints in 
this system are combinations of the derivatives of the 
dynamical fields. So we normalize the curves in this figure 
by the quantity | |Vu(i)| |, which is the natural coordinate- 
invariant L 2 measure of the derivatives of the dynamical 
fields: 

||V M (t)|| 2 = fg^V i u a V j u p S aP g L ^d a x. (71) 

The ratio of these quantities, ||C(t)||/||Vu(f)||, is there- 
fore a meaningful dimcnsionlcss measure of the magni- 
tude of constraint violations. When the value of this 
ratio becomes of order unity, the dynamical fields do not 
satisfy the constraints at all. As we can see in Fig. ^ 
the constraint satisfying initial data quickly evolve to a 
state in which this constraint measure is of order unity. A 
large increase in constraint violation occurs as the outgo- 
ing scalar wave pulse passes through the outer boundary 
of the computational domain. After this time the nu- 
merical solution to the first-order scalar wave system no 
longer represents a solution to the original scalar field 
equation. 

In Fig.|21we demonstrate that these numerical solutions 
are nevertheless numerically convergent. We measure the 
convergence of these solutions by depicting the quantity 

\\6u(t)\\ 2 = J S aP {u% r -u a R ) g l ' 2 d s x, 

(72) 
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FIG. 3: Constraint violations for evolutions with 71 = 72 = 0, 
constraint preserving boundary conditions, and no constraint 
projection. Solid curves are normalized by the quantity 
||Vw(t)|| while the dashed curves are normalized by ||Vu(0)||. 
Decay of the normalization factor ||Vw(t)|| rather than growth 
of the constraints causes the growth in the highest-resolution 
solid curves, which have constant roundoff- level constraint vi- 
olations. 
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FIG. 4: Convergence of evolutions shown in Fig. |^1 Plot- 
ted are differences from the evolution with N = 81, which is 
henceforth the reference solution iir. Solid curves are nor- 
malized by while the dashed curves are normalized 
by ||u(0)||. Decay of the normalization factor causes 
the growth in the highest-resolution solid curves, for which 
|[<kt(t)|| is constant at roundoff level. 



divided by a suitable normalization. This quantity mea- 
sures the difference between the solution obtained 
with radial resolution N r , compared to a reference solu- 
tion u^. In Fig. |2 we use the numerical solution com- 
puted with the largest number of radial basis functions 
(N r = 81 in this case) as the reference solution. In order 
to make these difference measures meaningful, we nor- 
malize them by dividing by an analogous measure of the 
solution itself: 

\\u(t)\\ 2 = f S a(3 u%y Nr gV 2 d 3 x. (73) 

Figure [21 shows that our computational methods are nu- 
merically convergent, even if the solutions are constraint 
violating and are therefore unphysical. The rate of con- 
vergence of these solutions changes at about t — 1QM be- 
cause a short wavelength reflected pulse enters the com- 
putational domain at about this time. The convergence 
of these solutions shows that these constraint violations 
are a feature of the evolution system and the boundary 
conditions, rather than being artifacts of a poor numeri- 
cal technique. 

Next we evolve the same initial data, Eqs. 
using the same standard scalar wave evolution equations 
(71 — 72 = 0) , but this time we use constraint preserving 
boundary conditions on the fields Z 1 and Zf, Eqs. if4"0"|) 
and (|41|l . We use no constraint projection in these evo- 
lutions. Figure |21 shows that the constraints are in fact 
satisfied by these solutions to truncation level errors. The 
solid curves in Fig. show the ratio ||C(t)||/||Vu(t)|| 



while the dashed curves show ||C(t)||/||Vu(0)||. The only 
difference is that the denominator used for the dashed 
curves is time independent. The solid curves show that 
the relative constraint error is approximately constant in 
time until about t = 40, at which time a truncation error 
level constraint-violating pulse from the outer boundary 
has advected inward across the grid and fallen into the 
black hole. After t = 40 the relative constraint error 
decreases with time. The highest-resolution solid curves 
behave differently: they increase exponentially with time. 
However, this growth occurs only because the normaliza- 
tion factor in the denominator (which measures the size 
of the derivatives of the fields) goes to zero as the scalar 
wave pulse leaves the computational domain. The high- 
est resolution dashed curves show that the absolute con- 
straint error for these resolutions is constant at roundoff 
level. 

Figure 01 illustrates the numerical convergence of these 
evolutions. Plotted are the ratios of the differences 
||<5w(t)|| to a measure of the size of the fields. The solid 
curves in Fig.0|show the ratio ||<$u(t)||/||u(t)|| while the 
dashed curves show ||<Jti(i)||/||it(0)||. Again, the only 
difference is that the denominator used for the dashed 
curves is time independent. Figures |31 and ^ show that 
these scalar field evolutions arc stable, constraint preserv- 
ing and numerically convergent. These solutions there- 
fore represent what we expect to be the correct physical 
solution to this problem. Were this our only objective, 
this paper would end here. However our primary interest 
here is to study the use of projection methods to control 
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the growth of constraints. So we will use the solution 
found here as a reference to which our later evolutions 
using constraint projection can be compared. 

Our last simulation to study the effects of boundary 
conditions on the growth of the constraints uses a non- 
standard scalar field evolution system with 71 = and 
72 = —l/M. In other respects, however, this simulation 
is identical to the one depicted in Figs. and 21 It uses 
the same initial data, Eqs. (|60|I - H62(I . the same constraint 
preserving boundary conditions, and no constraint pro- 
jection. Because we use Eq. I)39fl as a boundary condition 
on U l ~ , the constraint-preserving solutions of the equa- 
tions are the same as those obtained with 71 =72 = 0. 
However, using an evolution system with 72 = — l/M 
introduces unstable bulk terms into the constraint evo- 
lution equations, Eqs. Q29J) and (|30|l . so the constraint- 
violating solutions of the equations will be different. Con- 
sequently this system is much more pathological than 
the standard scalar field system, and provides a much 
more difficult challenge for the constraint control meth- 
ods studied here. Figure |3] shows the evolution of the 
constraints in this system. Truncation level constraint 
violations in the initial data grow exponentially with an 
c-folding time of approximately 1.1M in these evolutions. 
The ratio ||C(f)||/||Vu(£)|| approaches a constant of or- 
der unity at late times once the constraint violating por- 
tion of the solution dominates and the denominator be- 
gins to grow exponentially as well. The small inset graph 
in Fig.[S]illustrates that the divergence of these solutions 
from the reference solution of Fig. 0] grows at the same 
rate for all spatial resolutions. This suggests that the 
growth is caused by a constraint violating solution to the 
evolution equations rather than a numerical instability. 

These evolutions with 72 = — l/M demonstrate that 
constraint preserving boundary conditions alone are in- 
sufficient to control the growth of constraints in this 
system. Since the Einstein evolution system is also be- 
lieved to contain bulk generated constraint violations , 
this example suggests that constraint preserving bound- 
ary conditions alone will not be sufficient to control the 
growth of the constraints in the Einstein system. 



B. Testing Constraint Projection 

In this section we discuss two numerical evolutions that 
explore the use of the constraint projection methods de- 
veloped in Sees. UTI and 1111 CI The first evolution uses the 
standard scalar wave evolution system with 71 = 72 = 0, 
and freezing boundary conditions. We have already seen 
in Figs. ^ an d ID that such evolutions exhibit significant 
constraint violations once the scalar wave pulse passes 
through the outer boundary of the computational do- 
main. In this numerical experiment we freely evolve the 
scalar field to the time t = 20M, and then perform a sin- 
gle constraint projection on the solution using Eqs. I|52|l - 
(|54|l with A = 2/M. We then evolve the system freely 
again to t — 40M. Figure El shows how the constraints 
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FIG. 5: Constraint violations for evolutions with 72 = —l/M, 
constraint preserving boundary conditions, and without con- 
straint projection. The inset shows differences [[<Su(t)|[/|[«(£)[| 
from the reference solution of Fig. 0] The curves level off at 
late times because both numerator and denominator grow ex- 
ponentially at the same rates. 



respond to a single constraint projection. We use a very 
fine time scale in Fig. showing in detail the times 
around t = 20M when the constraint projection is per- 
formed. Individual points in Fig. show the amount of 
constraint violation after each individual time step. The 
value of the constraints drops sharply at the time step 
where the constraint projection is performed, and as we 
expect, the value of the constraints after this projection 
step is smaller for higher resolutions. So the constraint 
projection step is successful in significantly reducing the 
size of the constraints. But something rather unexpected 
happens next: the constraints increase by orders of mag- 
nitude on the very next free evolution time step after the 
constraint projection. The small inset in Fig.|Hlshows the 
same data plotted on a linear rather than a logarithmic 
scale. This shows that the constraints grow linearly in 
time after the constraint projection step on a very short 
time scale. 

Figure [3 provides some information about the reason 
for this strange behavior by showing the convergence of 
these numerical solutions. For times before the constraint 
projection step at t = 20M, the solutions show good nu- 
merical convergence as the number of radial collocation 
points is increased. But there is a sharp breakdown of nu- 
merical convergence (or at least a sharp drop in the rate 
of numerical convergence) after the constraint projection 
step. 

Figure |S1 provides some deeper insight into the reason 
for this lack of convergence. Plotted in Fig. [H] are a se- 
quence of curves showing the radial dependences of the 
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FIG. 6: Constraint violations for evolutions with 71 = 
72 = 0, freezing boundary conditions, and a single con- 
straint projection at t — 20M (with A = 2/M). Points show 
I |C(t)| l/l I Vu(t)| I after each time step. The inset plots the 
same data on a linear scale. 



FIG. 8: Radial profiles of (^>)io and {CiC l )oo for the evolution 
of Fig. El The solid lines represent times t/M = 20, ... , 25. 
The dashed line represents the state just before the constraint 
projection at t/M = 20. The arrows indicate the location of 
the non-smoothness in ip. 




FIG. 7: Convergence of evolutions shown in Fig. |S| Plotted 
are differences from the evolution with N r = 81. 



dipole part of the scalar field (ip) 10 and the monopole 
part of the constraints (CiC 1 ) 00 at a sequence of times 
including the constraint projection step. The spherical 
harmonic components of a function Q are defined by 



(Q)i m = J Yi^ n (e,ip)Q(r,0, ( p) sinO dddip. (74) 
The dashed lines at the bottom of Fig. [S] shows the radial 



profiles at t = 20M immediately before the constraint 
projection, while the lowest solid lines show these pro- 
files at the same time t = 20M just after the projection. 
We see that the constraints essentially vanish after the 
constraint projection step. The next profile at t = 21M 
shows that the scalar field develops some non-smooth 
radial structure immediately after the projection step, 
which subsequently propagates into the computational 
domain. This non-smoothness in ip causes a sharp spike 
in the constraints, seen clearly in Fig. |H1 Spectral meth- 
ods do not converge well for non-smooth functions, so the 
emergence of this structure in ip explains the breakdown 
in the numerical convergence and thence the breakdown 
in our constraint projection method. The emergence of 
the non-smoothness in ip seems to be caused by the con- 
straint projection step in the following way: The projec- 
tion produces a ip that is non- vanishing at the boundary, 
and the freezing boundary condition then forces tp = Z 1 
(and Zf) to develop kinks (see Ref.Q) which propagate 
into the computational domain during the free evolution 
steps following the projection. 

Figures EHE1 demonstrate that constraint projection is 
not successful in removing large constraint violations 
when used in conjunction with freezing boundary con- 
ditions. One might hope that this failure could be cor- 
rected by projecting out the constraints before they are 
allowed to grow too large. Figure shows the conver- 
gence of solutions in which a constraint projection is 
performed after each evolution time step, for a variety 
of different time steps At. Like the evolutions shown 
in Figures EHH1 these evolutions use the standard scalar 
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FIG. 9: Differences between evolutions with time step At and 
the reference solution ur (of Fig. at fixed evolution times 
to. Evolutions use 71 = 72 = 0, freezing boundary conditions, 
and constraint projection with A = 2/M after each time step, 
AT = At. 



field system (71 =72 = 0), freezing boundary conditions, 
and constraint projection with A = 2/M. The three 
curves in Fig. measure the convergence of the solution 
(relative to the highest resolution reference solution de- 
picted in Fig.|3J at three different times in this evolution, 
to = 10.24M, 20.48M, and 30.72M. All of these evolu- 
tions use the same spatial resolution, N r = 51. These 
graphs show that the convergence towards the reference 
solution is only first order in the time step At. This 
convergence is significantly worse than that expected for 
the fourth-order Rungc-Kutta time step integrator that 
we use. In contrast the free evolutions with constraint 
preserving boundary conditions shown in Fig. 0] achieve 
||#u(io)||/||ti(to)|| ^ 10~ 10 with a timestep similar to the 
largest At shown in Fig. [5] We conclude that constraint 
projection produces only first-order in time convergent 
numerical solutions when used in conjunction with stan- 
dard freezing boundary conditions, and is therefore an 
ineffective substitute for constraint preserving boundary 
conditions. 

Finally we apply constraint projection to the patho- 
logical scalar wave evolution system (71 = and 72 = 
— l/M), which we failed to control with constraint pre- 
serving boundary conditions alone. We project every 
AT = 2M using A = y/2/M, and we continue to use 
constraint preserving boundary conditions. Except for 
constraint projection, this is the same as the evolution 
shown in Fig. [5] Figure ^3 shows that the constraints 
are reduced to truncation error levels in these evolutions. 
The small inset graph shows these same curves with a 
finer time resolution, so the saw-tooth shaped evolution 
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FIG. 10: Constraint violations ||C(t)||/||Vu(t)|| for evolu- 
tions with 71 = and 72 = —l/M, constraint preserving 
boundary conditions, and constraint projection with A = 
y/2/M every AT = 2M. Inset shows the same data with 
finer time resolution. 
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FIG. 11: Differences from the reference solution ur (of Fig.|lJ 
for the evolutions shown in Fig. 1101 



of the constraints can be seen more clearly. We note that 
constraint projection does not occur at every evolution 
time step in these simulations, but rather at fixed times 
separated by AT = 2M. The evolutions with the finest 
spatial resolution take more than one thousand time steps 
between projections. Figure 1111 shows the convergence 
between these evolutions and the highest resolution ref- 
erence solution depicted in Fig. 0] This figure demon- 
strates that the constraint projection method combined 
with constraint preserving boundary conditions succeeds 
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FIG. 12: Differences from the reference solu- 

tion ur of Fig. [I] are plotted for different choices of A. Evo- 
lutions with 7i = and 72 = — 1/M, constraint preserving 
boundary conditions, constraint projection every AT = 2M. 




FIG. 13: Evolutions with 71 = and 72 = —1/M, constraint 
preserving boundary conditions and constraint projection ev- 
ery AT. Differences from the reference solution ur (of Fig.[lJ 
at to — lOOAf for different choices of AT and A. 



in producing the same numerical solution as our reference 
solution — even for this pathological scalar field system. 



C. Optimizing Constraint Projection 

In this section we explore ways to optimize the use of 
the constraint projection methods developed in Sees. |n] 
and IIII CI In particular we investigate how important 
the choice of the parameter A is to the effectiveness of 
the projection, and we determine its optimal value. We 
also vary the time between projection steps, AT, and 
determine the optimal rate at which to perform these 
projections. Finally we measure the computational cost 
of performing a scalar field evolution with constraint pro- 
jection, compared to the cost of doing a free evolution. 

Figurc^Jshows convergence plots for evolutions of the 
pathological scalar field system with 71 = and 72 = 
— 1/M, constraint preserving boundary conditions, and 
constraint projection every AT = 2M. All evolutions use 
the same radial resolution, N r = 41. Each of the solid 
curves in Fig. II 21 represents an evolution using a different 
choice of the parameter A. We see that the evolutions 
using projections with A = \/2/M are somewhat closer 
to the reference solution than the others, but the size of 
the differences are not very sensitive to the value of A. 
The only projected solution having significantly worse 
accuracy than the others is the one with A = 00, which 
corresponds to the simple projection with ip = ip, II = 
fl and <I>i = diip. For all choices of A, including A = 
00, these evolutions are exponentially convergent with 



increasing A^. 

We have some understanding of why there is an op- 
timal choice for the parameter A: It is possible to ana- 
lyze the projection process completely and analytically 
for scalar field evolutions with a flat background metric 
on a computational domain with three-torus (T 3 ) topol- 
ogy. By performing a Fourier transform of the fields in 
this case it is easy to show that the fields break up into 
modes that propagate with the usual dispersion relation 
lo 2 = k ■ k, plus others that grow exponentially in time 
with dispersion relation 10 = 172- The projection step be- 
comes a simple algebraic transformation on the Fourier 
components of the field in this case. So it is straightfor- 
ward to show that the projection step completely removes 
the modes that grow exponentially with time only when 
the parameters satisfy A 2 = 27! . For evolutions on com- 
putational domains with different topologies, and differ- 
ent background metrics, it is not possible to determine 
the optimal choice of A using such a simple argument. 
However it is not surprising that the optimal choice is 
not too different from A 2 = 2t| . 

Next we consider the effect of varying the times be- 
tween constraint projections. Figure ITS1 shows the con- 
vergence measure ||<$u(io)||/||u(£o)|| for evolutions of the 
pathological scalar field system with 71 = and 72 = 
— 1/M, constraint preserving boundary conditions, and 
constraint projections with various values of A and AT. 
These evolutions are all carried out with the same radial 
resolution N r = 41, and are compared with the reference 
solution of Fig. 0] at the time to = 100M. Each curve 
in Fig. E| represents a set of evolutions with the same 
value of A but varying AT. The smallest AT for each 
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FIG. 14: Solid curve (left axis) shows the ratio of time spent in 
elliptic solves to time spent in the hyperbolic evolution code. 
Dashed curve (right axis) shows the ratio of time required for 
one elliptic solve to the time for one evolution time step. 

curve corresponds to projecting at each evolution time 
step. We see that all of these curves show a minimum 
difference with the reference solution, and this minimum 
occurs at about AT « 1M in all of these curves. This 
coincides with the e-folding time of the bulk constraint 
violations, — 1/72; hence we expect that constraint pro- 
jection should generally be applied on a time-scale com- 
parable to that of the constraint growth. Figure IT31 also 
reveals that projections performed with A 2 = 27! are 
the optimal ones over a fairly broad range of projection 
times AT. The evolutions with simple constraint projec- 
tion (A = 00) crash for very small values of AT, as well 
as for AT = 10M. 

Finally, we have made some measurements to evalu- 
ate the computational cost of doing scalar field evolu- 
tions with constraint projection, compared to the cost of 
free evolution. Figure PHI contains two curves that mea- 
sure the computational cost of doing optimal projection 
with AT = 2M. The solid curve shows the ratio of the 
time the code spends doing the constraint projection step 
(i.e. doing the elliptic solve) T n with the time the code 
spends doing evolution steps Th yp . This ratio decreases 
from about 0.1 using a very coarse spatial resolution to 
about 0.003 using a very fine spatial resolution. The 
ratio T n/Thyp decreases when the spatial resolution is 
increased because the code must take many more free 
evolution time steps in the time AT between projection 
steps in this case. The dashed curve in Fig. 1141 measures 
the ratio of the time needed to perform one constraint 
projection, t e \\, with the time needed to take one free 
evolution step, ihyp- We see that this ratio is fairly in- 
dependent of resolution using our spectral elliptic solver, 
and ranges from about 3.5 at low spatial resolution to 



about 5 at high resolution. These tests show that the 
computational cost of performing constraint projection 
is only a small fraction of the total computational cost 
of performing these scalar field evolutions. We conclude 
that computational cost should not be used as an argu- 
ment against the use of constraint projection methods. 



V. DISCUSSION 

We have developed general methods in Sec. El for con- 
structing optimal projection operators that map the dy- 
namical fields of hyperbolic evolution systems onto the 
constraint submanifold associated with these systems. 
These methods are worked out explicitly in Sec. II I II for 
the case of a new evolution system that describes the 
propagation of a scalar field on a fixed background space- 
time. The constraint projection map for this system re- 
quires the solution of one elliptic partial differential equa- 
tion each time a projection is performed. The new scalar 
field system introduced in Sec. IIIII has the interesting 
property that it suffers from constraint violations that 
flow into the domain through timelike boundaries and 
also from violations generated by bulk terms in the equa- 
tions. So this system exhibits both types of constraint 
violating pathologies that can occur in the Einstein evo- 
lution system. To test our constraint projection meth- 
ods, we perform a number of numerical evolutions of this 
scalar field system propagating on a black-hole space- 
time. We show that constraint preserving boundary con- 
ditions alone are not capable of controlling the growth 
of constraints in this scalar field system. Constraint pro- 
jection is also shown to be ineffective when used in con- 
junction with traditional boundary conditions that do 
not prevent the influx of constraint violations through 
the boundary. However we show that the combination of 
constraint projection and constraint preserving bound- 
ary conditions is a very effective method for controlling 
the growth of the constraints. We measure the compu- 
tational cost of performing these constraint projections 
and show that at the highest numerical resolutions, the 
projections account for only a fraction of a percent of 
the total computational cost of the evolution. Thus high 
computational cost can no longer be cited as a reason to 
avoid constraint projection techniques. 
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